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Conversion  Factors 

SI  to  Inch/Pound 


Multiply 

By 

To  obtain 

Length 

centimeter  (cm) 

0.3937 

inch  (in.) 

millimeter  (mm) 

0.03937 

inch  (in.) 

meter  (m) 

3.281 

foot  (ft) 

Area 

square  centimeter  (cm1 2) 

0.001076 

square  foot  (ft2) 

Volume 

milliliter  (mL) 

0.0338 

ounce 

liter  (L) 

0.2642 

gallon  (gal) 

cubic  centimeter  (cm3) 

0.06102 

cubic  inch  (in3) 

Mass 

gram  (g) 

0.0022 

pound  (lb) 

Density 

gram  per  cubic  centimeter  (g/cm3) 

62.4220 

pound  per  cubic  foot  (lb/ft3) 

Hydraulic  conductivity 

centimeter  per  day  (m/d) 

0.3937 

inch  per  day  (in./d) 

Concentration 

milligram  per  liter  (mg/L) 

0.000062426 

pound  per  cubic  foot  (lb/ft3) 

microgram  per  liter  (pg/L) 

0.0062426 

pound  per  cubic  foot  (lb/ft3) 

milligram  per  kilogram  (mg/kg) 

0.002 

pound  per  ton  (Ib/ton) 

Temperature  in  degrees  Celsius  (°C)  may  be  converted  to  degrees  Fahrenheit  (°F)  as  follows: 

°F=(1.8x°C)+32 
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CA  -  California 
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ET  -  Evapotranspiration 

GLEAMS  -  Groundwater  Loading  Effects  of 
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GUI  -  graphical  user  interface 
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culture 
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Abstract 

Seven  unsaturated-zone  solute-transport  models  were 
tested  with  two  data  sets  to  select  models  for  use  by  the 
Agricultural  Chemical  Team  of  the  U.S.  Geological  Survey’s 
National  Water-Quality  Assessment  Program.  The  data  sets 
were  from  a  bromide  tracer  test  near  Merced,  California,  and 
an  atrazine  study  in  the  White  River  Basin,  Indiana.  In  this 
study  the  models  are  designated  either  as  complex  or  simple 
based  on  the  water  flux  algorithm.  The  complex  models, 
HYDRUS2D,  LEACHP,  RZWQM,  and  VS2DT,  use  Richards’ 
equation  to  simulate  water  flux  and  are  well  suited  to  process 
understanding.  The  simple  models,  CALF,  GLEAMS,  and 
PRZM,  use  a  tipping-bucket  algorithm  and  are  more  amenable 
to  extrapolation  because  they  require  fewer  input  parameters. 
The  purpose  of  this  report  is  not  to  endorse  a  particular  model, 
but  to  describe  useful  features,  potential  capabilities,  and  pos¬ 
sible  limitations  that  emerged  from  working  with  the  model 
input  data  sets.  More  rigorous  assessment  of  model  applicabil¬ 
ity  involves  proper  calibration,  which  was  beyond  the  scope  of 
this  study. 

Uncalibrated  (“cold”)  simulations  were  run  using  all 
seven  models  to  predict  the  transport  of  bromide  (Merced) 
and  the  transport  and  fate  of  atrazine  and  three  of  its  transfor¬ 
mation  products  (White  River  Basin).  Among  the  complex 
models,  HYDRUS2D  successfully  predicted  both  the  surface 
retention  and  accumulation  of  bromide  at  depth  at  the  Merced 
site,  whereas  RZWQM  and  VS2DT  predicted  only  the  latter. 
RZWQM  predictions  of  atrazine  were  closest  to  observed  val¬ 
ues  at  the  White  River  Basin  site,  where  preferential  flow  has 
been  observed.  LEACHP  predicted  smaller  solute  concentra¬ 
tions  than  observed  at  both  the  Merced  and  White  River  Basin 
sites.  Among  the  simple  models,  CALF  predicted  the  highest 
values  of  atrazine  and  deethylatrazine  at  the  measurement 
depth  of  1.5  meters.  CALF  includes  the  Addiscott  flow  option 
for  preferential  flow,  and  also  accepts  user-specified  disper- 
sivity.  PRZM  underpredicted  solute  concentrations,  probably 
because  control  of  dispersion  is  a  problem  with  this  model. 
GLEAMS  has  a  maximum  simulation  depth  of  1.5  meters, 
which  is  limiting  for  mass-balance  purposes  because  it  creates 
a  potential  disconnect  between  unsaturated-zone  transport  and 
the  water  table. 


Of  the  models  tested,  RZWQM,  HYDRUS2D,  VS2DT, 
GLEAMS  and  PRZM  had  graphical  user  interfaces.  Extensive 
documentation  was  available  for  RZWQM,  HYDRUS2D,  and 
VS2DT.  RZWQM  can  explicitly  simulate  water  and  solute 
flux  in  macropores,  and  both  HYDRUS2D  and  VS2DT  can 
simulate  water  and  solute  flux  in  two  dimensions.  The  version 
of  RZWQM  tested  had  a  maximum  simulation  depth  of  3 
meters.  The  complex  models  simulate  the  formation,  transport, 
and  fate  of  degradates  of  up  to  three  to  five  compounds  includ¬ 
ing  the  parent,  with  the  exception  of  VS2DT,  which  simulates 
the  transport  and  fate  of  a  single  compound. 


Introduction 

One  of  the  primary  goals  in  the  second  decade  of  the  U.S. 
Geological  Survey’s  (USGS)  National  Water-Quality  Assess¬ 
ment  Program  (NAWQA)  is  to  refine  understanding  of  key 
processes  that  control  contaminant  behavior  in  the  Nation’s 
water  resources.  In  particular,  extensive  use  of  fertilizers  and 
pesticides  in  agricultural  areas  can  affect  the  quality  of  surface 
and  ground  waters.  Improved  understanding  of  relations 
among  contaminant  sources,  their  transport  through  the  hydro- 
logic  system,  and  effects  on  water  quality  will  provide  the 
basis  for  predicting  water-quality  conditions  in  unmonitored 
areas.  To  achieve  this  goal,  the  Agricultural  Chemical  Team 
(ACT)  is  conducting  studies  of  small  stream  basins  (tens  of 
square  kilometers  in  size)  in  five  agricultural  areas  associated 
with  selected  2001  NAWQA  study  units.  The  ACT  studies 
focus  on  the  major  aspects  of  the  hydrologic  cycle,  including 
the  unsaturated  zone,  ground-water  flow,  and  ground-water/ 
surface-water  interactions. 

Although  numerous  unsaturated-zone  (UZ)  model 
comparison  studies  have  been  conducted  by  other  research¬ 
ers  in  specific  locales  (appendix),  the  NAWQA  ACT  stud¬ 
ies  encompass  a  diverse  range  of  environmental  settings 
and  management  practices  across  both  semiarid  and  humid 
climates.  Because  of  this  diversity,  both  one-dimensional 
and  two-dimensional  flow  regimes  are  expected,  and  some 
preferential  flow  is  anticipated.  Currently,  no  single  UZ  code 
accommodates  all  of  these  conditions.  Several  models  were 
tested  with  available  data  to  determine  potential  applicabil- 
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ity  to  the  NAWQA  data  sets.  The  four  “complex  models” 
that  were  tested  use  Richards’  equation  for  water  flux.  These 
models  will  be  used  to  evaluate  physical  and  chemical  pro¬ 
cesses  affecting  contaminant  transport  and  fate  at  intensive 
UZ  monitoring  sites.  Complex  models  typically  require  more 
input  data  than  their  simple  counterparts.  In  contrast,  the  three 
“simple  models”  tested  use  a  “tipping-bucket”  algorithm  for 
water  flux,  and  therefore,  require  fewer  input  data.  The  tipping 
bucket  algorithm  assumes  that  water  in  excess  of  field  capacity 
moves  down  to  the  next  soil  compartment,  provided  that  the 
next  compartment  is  not  saturated.  Simple  models  will  be  used 
to  extrapolate  data  from  intensive  UZ  study  sites  to  locations 
with  limited  data  (e.g.,  soil  texture  from  soil  survey  databases). 
The  simple  model  results  will  be  used  to  derive  mass  balances 
of  water  and  solutes  for  the  entire  stream  basin  in  each  ACT 
study  area. 

Purpose  and  Scope 

The  objective  of  this  study  was  to  evaluate  ease  of  use 
and  differences  in  model  output  as  related  to  model  capabili¬ 
ties.  Differences  in  model  output  have  potential  implications 
for  model  applicability,  but  are  difficult  to  evaluate  without 
calibrating  the  models.  Calibration  can  reduce  model  param¬ 
eter  uncertainty  and  allow  for  more  robust  comparison  of  mod¬ 
eling  results.  Model  calibration  was  not  performed,  however, 
because  this  was  outside  of  the  scope  of  this  study.  Even  with¬ 
out  calibration,  differences  in  model  output  can  yield  insight 
into  each  model’s  ability  to  represent  the  physical  system  and 
processes  affecting  that  system.  For  example,  certain  tipping- 
bucket  models  cannot  simulate  upward  water  flux  by  capillary 
rise  so  they  tend  to  overpredict  percolation  (Garratt  and  others, 
2002).  Similarly,  models  that  do  not  simulate  preferential 
flow  might  be  inappropriate  in  areas  with  highly  structured, 
cracked  soils  (such  as  the  White  River  Basin,  Indiana),  and 
their  predictions  often  differ  from  dual  porosity  models  such 
as  RZWQM. 

The  purpose  of  this  report  is  not  to  endorse  a  particular 
model,  but  to  describe  the  results  of  the  model  evaluations. 
Features  of  20  UZ  models  were  compiled  in  a  spreadsheet  for 
initial  comparison.  Seven  were  selected  for  evaluation  using 
test  data  sets  (table  1),  based  on  the  following  criteria: 

•  Quality  of  documentation 

•  Ease  of  use 

•  Available  within  the  public  domain 

•  Source  code  available  for  modification 

•  History  of  successful  use  by  regulatory  or  other  agency 

•  Anticipation  of  ongoing  support 

•  At  least  one  complex  model  for  process  understanding 

•  One  simple  model  for  scale-up  and  mass-balance  esti¬ 
mates 


•  At  least  one  complex  model  should  have  preferential 
flow  capability 

•  At  least  one  complex  model  should  have  two-dimen¬ 
sional  (2-D)  capability 

The  complex  models  chosen  for  evaluation  were 
HYDRUS2D,  FEACHP,  RZWQM,  and  VS2DT;  and  the 
simple  models  were  CAFF,  GFEAMS,  and  PRZM.  Major 
processes  simulated  by  the  models  are  shown  in  table  1 . 

The  test  data  sets  for  model  evaluation  were  selected 
using  several  criteria  established  on  the  basis  of  anticipated 
UZ  modeling  needs  in  ACT  study  areas.  These  criteria  were: 

•  Comprehensive  field  data  set  for  model  input  and  com¬ 
parison  with  model  output,  including  on-site  climate 
data,  soil  descriptions,  and  soil-moisture  and  solute- 
concentration  measurements  through  time  and/or 
space. 

•  Both  conservative  transport  (e.g.,  bromide)  and  reactive 
transport  with  parent  compound  (e.g.,  atrazine)  and 
degradates  in  measurable  concentrations. 

•  Both  relatively  homogeneous  soils  and  also  structured 
soils  with  documented  preferential  flow. 

•  Range  of  climate  conditions. 

•  Focated  near  ACT  studies  of  interest. 

The  two  data  sets  satisfying  these  criteria  that  were  used 
in  the  model  evaluation  were: 

•  Merced,  California  (CA),  bromide  tracer  concentration 
profile  measured  January  25,  1993  (located  near  the 
San  Joaquin-Tulare  River  Basins  ACT  study);  and 

•  White  River  Basin,  Indiana  (IN),  atrazine  study,  site 
C,  sample  depth  150  centimeters  (cm)  (data  collected 
May  -  September  1994)  (located  near  the  White, 

Great,  and  Fittle  Miami  River  Basins  ACT  study). 

The  Merced,  CA,  site  features  a  conservative  tracer 
(bromide)  and  the  White  River  Basin,  IN,  site  features  reactive 
compounds,  including  measurable  concentrations  of  degra¬ 
dates.  The  Merced  soils  are  relatively  homogeneous  whereas 
the  White  River  Basin  site  has  highly  structured  soils  exhibit¬ 
ing  preferential  flow  (Bayless,  2001).  The  sites  encompass 
both  semiarid  (CA)  and  humid  (IN)  climate  conditions. 

Because  of  time  and  budget  constraints,  modeling  was 
limited  to  one  lysimeter  location  and  depth  at  the  White  River 
Basin  site,  and  one  concentration  profile  on  a  single  date  in  the 
case  of  the  Merced  site.  Five  modelers  simultaneously  tested 
one  or  more  of  the  seven  models  so  that  the  evaluations  could 
be  conducted  in  a  timely  manner.  The  bulk  of  the  time  spent 
on  the  project  involved  becoming  familiar  with  the  codes  and 
parameterizing  the  models  for  the  two  sites.  Concerted  effort 
was  dedicated  to  acquiring  an  internally  consistent  set  of  pes¬ 
ticide  parameters  that  reflected  the  literature  and  emphasized 
the  most  recently  published  information. 
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4  Evaluation  of  Unsaturated-Zone  Solute-Transport  Models  for  Studies  of  Agricultural  Chemicals 


Approach 

Each  of  the  seven  unsaturated  zone  models  was  tested 
with  the  two  data  sets  and  evaluated  for  ease  of  use  and  poten¬ 
tial  to  simulate  processes  of  interest.  All  models  were  tested 
via  “cold”  or  forward  simulation,  without  calibration.  Proper 
calibration  involves  inverse  modeling  to  estimate  parameter 
values  that  optimally  fit  the  observed  data  (Hill  and  others, 
2000).  Manual  calibration  was  not  considered  because  the 
process  entails  numerous  disadvantages,  such  as  subjectivity 
in  visual  assessment  of  fit  and  in  deciding  when  to  end  the 
calibration;  the  difficulty  of  simultaneously  calibrating  more 
than  two  parameters;  and  the  lack  of  statistical  information  on 
calibrated  parameters  (Dubus  and  others,  2002).  In  contrast, 
cold  simulation  involves  a  single  run  using  the  best  estimates 
of  model  input  parameters.  Cold  simulation  tests  both  the 
model  and  the  parameterization  process,  which  frequently 
involves  informed  judgment  and  use  of  algorithms  to  derive 
parameters  not  readily  available  (e.g.,  water-retention  func¬ 
tions  in  complex  models).  Finally,  cold  simulation  is  appropri¬ 
ate  where  potential  use  without  calibration  is  anticipated,  such 
as  by  a  regulatory  agency. 

The  physical  domain  was  set  up  using  the  simulation 
layers  shown  in  table  2  and  corresponding,  average  values  of 
sediment  texture.  Depth  intervals  with  similar  texture  were 
merged  into  thick  layers  at  the  Merced  site.  At  the  White  River 
Basin  site,  all  layers  were  simulated  to  test  model  capabili¬ 
ties  in  heterogeneous,  highly  structured  soils.  Key  sediment 
properties  are  shown  for  each  UZ  sampling  interval  in  table 
2.  Saturated  hydraulic  conductivity  (K  j  was  obtained  from 


field  or  laboratory  measurements  at  both  sites.  Climate  data 
are  from  a  combination  of  onsite  data  and  data  from  nearby 
National  Oceanic  and  Atmospheric  Administration  (NOAA) 
stations. 

Objective  comparison  of  the  models  required  use  of  the 
same  input  values  for  parameters  common  to  all  seven  models. 
Compounds  simulated  by  the  models  were: 

•  Merced,  CA: 
bromide  (Br) 

•  White  River  Basin,  IN: 
atrazine  (ATR)  (parent) 
desethylatrazine  (DEA)  (daughter) 
hydroxyatrazine  (HYA)  (daughter) 

didealkylatrazine  (DDA)  (granddaughter) 

Values  of  fate  and  transport  parameters  common  to  the 
models  are  shown  in  table  3.  In  the  case  of  atrazine  and  its 
degradates,  best  estimates  based  on  literature  sources  generally 
were  used;  however,  transformation  rates  based  on  a  prior  field 
study  (Bayless,  2001)  were  used  in  lieu  of  laboratory  values 
from  the  literature.  The  latter  were  2-4  orders  of  magnitude 
lower  than  the  field-derived  values  (0.025  -  0.035  day1), 
which  would  have  significantly  reduced  predicted  concentra¬ 
tions  of  atrazine  degradates.  Field  transformation  rates  were 
emphasized  because  laboratory  conditions  do  not  represent 
transient  and  spatially  variable  conditions  encountered  in 
the  field.  Laboratory  transformation  rates  commonly  differ 


Table  2.  Simulation  layers  and  sediment  properties. 

[Kjat,  saturated  hydraulic  conductivity;  cm/d,  centimeters  per  day;  g/cm3,  grams  per  cubic  centimeter] 


Depth 

(cm) 

Simulated 

layer 

K  , 

sat 

(cm/d) 

Bulk  den¬ 
sity  (g/cm3) 

Sand 

Percent 

Silt 

Clav 

Merced.  CA 

5 

1 

37.4 

1.25 

76.29 

21.35 

2.36 

18 

1 

37.4 

1.68 

73.03 

24.43 

2.54 

37 

1 

37.4 

1.81 

74.39 

23.42 

2.20 

62 

1 

37.4 

1.83 

74.85 

23.13 

2.02 

92 

1 

37.4 

1.72 

75.86 

22.36 

1.78 

123 

2 

29.4 

1.77 

63.61 

33.16 

3.23 

151 

2 

29.4 

1.73 

55.15 

38.84 

6.01 

180 

3 

37.4 

1.83 

69.96 

27.18 

2.86 

212 

3 

37.4 

1.83 

88.46 

10.68 

0.85 

241 

3 

37.4 

1.82 

94.58 

4.98 

0.44 

257 

3 

37.4 

1.81 

93.71 

6.28 

0.01 

White  River  Basin,  IN 

15 

1 

35 

1.77a 

23.40 

49.50 

27.10 

28 

2 

35 

1.77 * 

14.80 

51.00 

34.20 

41 

3 

35 

1 .77“ 

13.90 

52.80 

33.30 

53 

4 

35 

1.77* 

16.80 

50.50 

30.20 

66 

5 

35 

1.77* 

29.00 

40.80 

30.20 

81 

6 

35 

1.77* 

32.10 

39.30 

28.60 

97 

7 

35 

1.77* 

31.90 

45.00 

23.10 

114 

8 

35 

1.77* 

33.60 

46.40 

20.00 

132 

9 

35 

1.77* 

36.10 

43.90 

20.00 

145 

10 

35 

1.77* 

39.60 

42.90 

17.50 

157 

11 

35 

1.77* 

34.80 

45.50 

19.70 

‘Average  at  site. 
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Table  3.  Model  parameters  associated  with  chemical  transport  of  bromide,  atrazine,  and  atrazine  degradates. 

[CA,  California;  IN,  Indiana;  ppm,  parts  per  million;  %,  percent;  mm1,  per  millimeter;  ml/g,  milliliters  per  gram;  Koc,  organic  carbon  partition  coefficient;  Kow, 
octanol-water  partition  coefficient;  mm  Hg,  millimeters  mercury;  cm2/sec,  square  centimeters  per  second;  cm,  centimeters;  NA,  not  applicable;  NU,  not  used 
because  atrazine  applied  before  plant  germination] 


Parameter 

Merced,  CA 

White  River  Basin,  IN 

Bromide 

Atrazine 

Desethylatrazine 

Hydroxyatrazine 

Didealkylatra- 

zine 

Water  solubility  (ppm) 

NA 

33 

2,700 

16 

94 

Foliar  residue  half-life  (days) 

NA 

NU 

NU 

NU 

NU 

Foliar  washoff  fraction  (%) 

NA 

NU 

NU 

NU 

NU 

Foliar  washoff  power  term 
(mm'1) 

NA 

NU 

NU 

NU 

NU 

Partition  coefficient,  Kd  (ml/g) 

NA 

1.75 

0.56 

4.26 

0.007 

Partition  coefficient  normal¬ 
ized  to  organic  carbon,  K 
(ml/g) 

NA 

250 

80 

609 

1 

Kow  (for  plant  uptake) 

NA 

370 

110 

950 

1 

Lumped,  overall: 

Transformation  rate  (day1)  in 
sorbed  phase 

NA 

0.085 

ATR  to  DEA:  0.025 
ATR  to  DIA:  0.035a 
ATR  to  HYA:  0.025 

DEA  to  DDA:  0.100 

Lumped,  overall: 

0.57 

Lumped,  over¬ 
all:  2.50 

Freundlich  isotherm  constant, 

Kf  (mL/g)b 

NA 

0.90 

0.57 

6.6 

0.97 

Freundlich  exponent  1/n 

NA 

0.85 

0.95 

0.76 

0.92 

Vapor  pressure  (mm  Hg) 

NA 

5.64E-06 

1.28E-04 

4.30E-05 

2.62E-02 

Henry’s  constant  (dimension¬ 
less) 

NA 

7.57E-08 

7.57E-08 

7.57E-08 

7.57E-08 

Diffusion  in  air  (cm2/ sec) 

NA 

0.058 

0.065 

0.060 

0.083 

Diffusion  in  water  (cm2/sec) 

5.79E-06 

4.4E-06 

5.0E-06 

4.6E-06 

6.3E-06 

Dispersivity  (cm) 

5 

360 

360 

360 

360 

aNot  simulated  in  this  study  because  observed  concentration  was  negligible,  but  contributes  to  lumped  atrazine  degradation  rate. 
bKf=KfocFoc 


from  field  rates  because  in  the  laboratory,  soil  samples  are 
disaggregated,  pesticide  solutions  are  uniformly  mixed,  and 
flow  rates  are  controlled  (Ma  and  Selim,  1996).  Additionally, 
laboratory  transport  studies  rarely  include  secondary  porosity 
as  a  dominant  feature,  or  the  non-steady  interactions  of  system 
chemistry  and  biology. 

Because  the  simulated  processes  varied  from  model  to 
model  (table  1),  each  has  certain  parameters  not  found  in  the 
other  models.  For  example,  RZWQM  requires  inputs  for  the 
dimensions  of  soil  cracks  and  the  “average  volume  fraction  of 
macroporosity.”  At  the  White  River  Basin  site,  these  values 
were  estimated  based  on  soil  cracks  observed  at  the  ground 
surface.  The  macropore  option  of  RZWQM  was  not  used  for 
the  Merced  site. 

The  complex  models  used  a  variety  of  characteristic 
curve  functions  to  specify  the  relations  between  moisture 


content,  matric  pressure,  and  hydraulic  conductivity.  The 
van  Genuchten  (1980)  function  was  used  for  HYDRUS  and 
VS2DT.  For  Merced,  VS2DT  simulations  used  a  modified 
version  of  the  van  Genuchten  function  that  included  sepa¬ 
rate  curve-fitting  parameters  for  retentivity  and  conductiv¬ 
ity  (Flealy  and  Ronan,  1996).  LEACHP  and  RZWQM  used 
Campbell’s  (Hutson,  2003)  and  Brooks-Corey  functions 
(Ahuja  and  others,  2000),  respectively. 

Parameters  were  derived  for  every  combination  of  sedi¬ 
ment  layer  (table  2)  and  characteristic  curve  function  for  both 
sites.  For  Merced,  the  van  Genuchten  retention  function  was 
fitted  to  field-measured  moisture  content  and  matric  pressure, 
and  the  conductivity  function  was  fitted  to  laboratory-mea¬ 
sured  moisture  content  and  unsaturated  hydraulic  conductivity. 
For  the  White  River  Basin,  grain-size  data  were  used  with  the 
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U.S.  Department  of  Agriculture’s  (USDA’s)  Rosetta  model 
(Schaap  and  others,  2001)  to  estimate  parameters  of  the  1980 
van  Genuchten  function  for  each  sediment  layer.  For  both 
sites,  parameters  of  Campbell’s  and  Brooks-Corey  functions 
were  optimized  using  the  Solver  tool  in  Microsoft  Excel  to 
achieve  characteristic  curves  as  similar  to  the  van  Genuchten 
functions  as  possible.  Figure  la  shows  the  fit  of  the  van 
Genuchten  function  for  sediment  type  1  at  Merced,  CA.  An 
example  of  the  Brooks-Corey  function  fitted  to  van  Genuchten 
values  from  the  White  River  Basin  ( 1 .45  meters  (m))  is  shown 
in  figure  lb.  To  construct  the  plots,  the  sign  convention  of  mat- 
ric  pressure  was  reversed  and  the  resulting  positive  values  (in 
centimeters)  were  log  transformed.  Moisture  contents  at  field 
capacity  (330  cm)  and  wilting  point  (15,000  cm)  were  taken 
from  the  van  Genuchten  curves  to  parameterize  the  simple 
models. 

Initial  conditions  of  pressure  and  solute  concentration 
were  specified  for  each  of  the  two  test  sites.  Initial  pressure 
was  set  for  Merced  based  on  field  data  using  a  linear  function 
that  ranged  from  -800  cm  at  0  cm  depth  to  -200  cm  at  600  cm 
depth.  The  initial  pressure  for  White  River  Basin  was  set  to  a 
constant  value  of  -10  cm.  For  White  River  Basin,  a  1 -month 
period  prior  to  the  application  of  pesticides  was  simulated 
to  establish  initial  conditions  for  solute  transport.  All  solute 
concentrations  were  initially  set  to  zero  at  both  sites. 

Top  boundary  conditions  consisted  of  weather  and  chemi¬ 
cal  inputs  to  the  different  models.  In  general,  the  upper  bound¬ 
ary  alternated  between  specified  flux  (infiltration)  and  model- 
computed  evapotranspiration.  Precipitation  records  from  the 
closest  daily-measured  raingage  were  used  to  calculate  the 
total  rainfall  and  timing  of  rain  events.  In  VS2DT,  infiltration 
was  applied  at  a  constant  rate  equal  to  the  mean  rainfall  rate 
for  the  study  period.  The  duration  of  infiltration  was  adjusted 
to  apply  the  correct  amount  of  precipitation  for  each  rainfall 
event.  In  the  remaining  models,  daily  rainfall  was  entered.  In 


Log  matric  pressure,  in  centimeters  (cm) 

Figure  la.  Fit  of  van  Genuchten  moisture-retention  function  to 
field  data  at  Merced,  California,  for  sediment  type  1  (0-125  cm  and 
175-300  cm  depths).  Parameters  of  the  van  Genuchten  function 
for  sediment  1  are:  0r  =  0.03;  tp  =  0.322;  |32  =  2.813;  a  =  0.017;  where 
0r  =  residual  moisture  content,  cp  =  porosity,  and  a  and  P2  are 
water-retention  parameters. 


general,  daily  precipitation  and  irrigation  records  were  used 
to  calculate  the  total  infiltration  while  solar  radiation,  wind 
movement  and  dew  point  data  were  required  to  calculate 
evapotranspiration . 

Bottom  boundary  conditions  were  specified  according 
to  the  capabilities  of  the  models.  In  particular,  LEACHP  was 
stable  with  White  River  Basin  data  only  when  water-table 
depth  was  fixed  at  the  bottom  of  the  soil  profile.  For  all  of 
the  models,  constant  matric  potential  was  specified  and  set  at 
-280  cm  for  Merced  (6  m  depth)  and  0  cm  for  the  White  River 
Basin  (2.5  m  depth). 


Results 

Differences  between  observed  and  simulated  values  can 
result  from  inappropriate  parameter  values  or  lack  of  the  mod¬ 
el’s  ability  to  simulate  processes  occurring  at  a  site.  Although 
the  models  in  this  study  were  not  calibrated  to  field  data,  dif¬ 
ferences  in  model  output  can  shed  light  on  model  capabilities. 
The  following  observations  are  based  on  visual  assessment  of 
model  fit  compared  with  observed  values  of  moisture  content, 
matric  potential,  and/or  chemical  concentrations.  Moisture 
content,  matric  potential,  and  bromide  concentration  were 
measured  at  the  Merced  site;  matric  potential  and  concentra¬ 
tions  of  atrazine  and  degradates  were  measured  for  the  White 
River  Basin  site. 


Merced,  CA 

Moisture-content  profiles  generated  by  complex  mod¬ 
els  generally  agreed  with  observed  values  at  depths  to  about 
120  cm  (fig.  2).  Below  this  depth,  LEACHP  and  RZWQM 


Log  matric  pressure,  in  centimeters  (cm) 

Figure  1b.  Fit  of  Brooks-Corey  moisture-retention  function  to  van 
Genuchten  values  derived  from  grain  size  data  at  the  White  River 
Basin  (1.45  meters).  Parameters  of  the  Brooks-Corey  function 
for  this  sediment  layer  are:  0r  =  0.043;  0s  =  0.316;  =  0.000614 

cm;  X  =  0.262;  B  =  0.672;  where  0r  =  residual  moisture  content,  0s 
=  saturated  moisture  content,  and  Aj;  X,  and  B  are  curve-fitting 
parameters. 
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Moisture  content,  in  cubic  centimeters  per 
cubic  centimeter  (cm3/cm3) 


Figure  2.  Moisture  content  with  depth  at  Merced,  California. 


more  closely  approximated  the  observed  profiles,  whereas 
HYDRUS2D  and  VS2DT  predicted  drier  conditions  than 
observed.  CALF  and  PRZM  generally  predicted  lower  mois¬ 
ture  content  than  observed,  which  may  have  resulted  from 
moisture  contents  at  field  capacity  (0.04  -  0.06)  and  wilt¬ 
ing  point  (0.03)  used  as  input  to  the  tipping-bucket  models. 
Although  these  values  seem  low,  they  correspond  to  an  empiri¬ 
cally  derived  water-retention  curve  for  the  Merced  site. 

GLEAMS  predictions  stopped  at  150  cm  because  that 
is  the  maximum  simulation  depth  for  this  model.  Also, 
GLEAMS  does  not  estimate  soil  moisture  per  se,  but  does  pro¬ 
vide  “soil  water  depth  by  layer”  in  centimeters.  These  values 
were  converted  to  tons  of  water  and  divided  by  the  tons  of  soil 
in  each  layer  to  produce  the  GLEAMS  output  in  fig.  2. 

HYDRUS2D,  LEACHP,  RZWQM,  and  VS2DT  ade¬ 
quately  simulated  matric  potential  to  about  125  cm  depth  (fig. 
3);  CALF  predicted  lower  potentials  than  observed  at  shal¬ 
low  depths.  HYDRUS2D,  RZWQM,  and  VS2DT  predicted 


decreasing  (i.e.,  more  negative)  matric  potential  below  125 
cm,  which  is  consistent  with  the  observed  values.  The  low 
pressures  at  depth  are  below  the  wetting  front  and  result  from 
drainage  and  evapotranspiration  (ET)  during  the  dry  spring 
and  summer  months  preceding  the  tracer  test.  LEACHP,  a 
complex  model,  predicted  nearly  constant  matric  potential 
(about  -100  cm)  with  depth. 

Bromide  tracer  was  applied  to  the  soil  surface  at  Merced 
on  August  20,  1992.  On  January  25,  1993,  observed  bro¬ 
mide  concentrations  at  Merced  were  bimodal,  with  some  of 
the  mass  retained  at  the  soil  surface  (less  than  50  cm)  (fig. 

4).  A  deeper  observed  peak  exists  at  about  160  cm.  Using  a 
mobile-immobile  feature,  HYDRUS2D  simulated  both  peaks, 
although  the  predicted  peak  concentration  of  mobile  bromide 
is  shallower  than  observed.  The  mobile-immobile  feature 
simulates  non-equilibrium  solute  transport  by  separating  the 
liquid  phase  into  mobile  (flowing)  and  immobile  (stagnant) 
zones,  and  solute  exchange  between  the  two  zones  is  modeled 
as  a  first-order  process  (Simunek  and  others,  1999).  RZWQM 
predicted  a  single  peak  at  depth  that  lagged  slightly  behind  the 
observed  values.  LEACHP  predicted  that  the  peak  bromide 
concentration  occurred  near  the  soil  surface.  Simple  models 
generally  underpredicted  bromide  concentration  at  depth. 
PRZM  predicted  the  lowest  bromide  concentrations  (less  than 
1  milligram  per  kilogram  (mg/kg));  concentrations  increase 
below  about  300  cm,  indicating  that  most  of  the  simulated 
bromide  peak  has  already  passed  through  the  model  profile. 
GLEAMS  simulated  a  concentration  peak  of  approximately  3 
mg/kg  at  the  1 00  cm  depth. 


White  River  Basin,  IN 

HYDRUS2D,  VS2DT,  RZWQM,  and  LEACHP  simu¬ 
lated  matric  potential  in  the  White  River  Basin  fairly  accu¬ 
rately  (fig.  5).  RZWQM  predictions  of  matric  potential 
decreased  (became  more  negative)  after  day  250,  but  there 
were  no  observed  values  at  later  times  for  comparison.  CALF 


EXPLANATION 
■  Measured 

-  VS2DT 

-  HYDRUS 

•  -  •  CALF 

- RZWQM 

-  LEACHP 


EXPLANATION 
■  Measured 

(avg.) 

-  VS2DT 

-  HYDRUS 

.  PRZM 

"  '  •  ‘  GLEAMS 

- RZWQM 

- LEACHP 

CALF 


Matric  potential,  in  centimeters-water 


Bromide  concentration,  in  milligrams  per  kilogram  (mg/kg) 


Figure  3.  Matric  potential  with  depth  at  Merced,  California. 


Figure  4.  Bromide  concentration  with  depth  at  Merced,  California. 
Units  are  milligrams  per  kilogram  of  soil  dry  weight. 
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EXPLANATION 
■  Measured 
matric  potential 

-  VS2DT 

-  HYDRUS 

- RZWQM 

- LEACHP 


Julian  days,  1994 

Figure  5.  Matric  potential  with  time  at  White  River  Basin,  Indiana. 

predictions  were  consistently  lower  (more  negative)  than 
observed  values.  PRZM  and  GLEAMS  do  not  simulate  matric 
potential. 

Atrazine  was  applied  on  April  18,  1994  (Julian  day  118), 
at  monitoring  site  C.  Bayless  (2001)  previously  argued  on  the 
basis  of  field  observations  and  modeling  results  that  preferen¬ 
tial  flow  is  important  at  this  site.  Use  of  the  macropore  option 
in  RZWQM  with  no  dispersivity  predicted  the  rapid  arrival  of 
atrazine  at  the  150-cm  sampling  depth  (fig.  6).  RZWQM  also 
simulated  atrazine  concentration  fairly  well  at  later  times  (150 
days  or  more).  Although  they  do  not  simulate  macropore  flow, 
HYDRUS2D  accurately  simulated  the  rapid  arrival  of  atrazine, 
and  VS2DT  predicted  atrazine  arrival  just  a  few  days  ahead  of 
observed  values.  The  rapid  transport  probably  resulted  from 
the  comparatively  large  dispersivity  value  used  for  these  simu¬ 
lations  (360  cm).  High  dispersivity  values  cause  the  leading 
edge  of  a  simulated  solute  plume  to  arrive  earlier,  increasing 
the  simulated  concentration  at  the  observation  point.  The  high 
dispersivity  value  was  taken  from  a  prior  modeling  study  in 
the  area  (Bayless,  2001).  Atrazine  concentrations  predicted 


by  HYDRUS2D  are  two  to  three  times  observed  values  after 
150  days.  Simulation  results  from  GLEAMS  and  PRZM  are 
presented  in  figures  6-8  for  a  depth  of  40  cm,  rather  than  1.5 
m,  because  both  models  predicted  negligible  concentrations 
of  atrazine  and  degradates  at  1.5  m  throughout  the  simulation 
period. 

HYDRUS2D  simulated  DEA  concentrations  fairly  accu¬ 
rately  up  to  about  160  days,  after  which  predicted  concentra¬ 
tions  were  about  twice  the  observed  values  (fig.  7).  RZWQM 
simulated  DEA  concentrations  fairly  accurately  after  about 
200  days.  All  other  models  (excluding  VS2DT,  which  does 
not  simulate  conversion  to  degradates)  underpredicted  DEA 
concentrations.  PRZM  predicted  extremely  low  DEA  concen¬ 
trations  at  40  cm  (10 7  micrograms  per  liter  (pg/L)  or  less), 
which  might  be  related  to  a  problem  with  the  grid  cell  size. 

The  grid  spacing  was  1 -millimeter  (mm)  layers  up  to  10  cm, 
and  1-cm  layers  at  depths  greater  than  10  cm,  as  suggested  by 
the  PRZM  documentation.  Grid  cell  size  influences  simulated 
numerical  dispersion  in  PRZM,  such  that  the  smaller  the  cell, 
the  lower  the  dispersion,  especially  for  a  strongly  sorbing  and 
decaying  solute.  Alternatively,  increasing  numerical  dispersion 
by  increasing  grid  size  would  accelerate  the  movement  of  the 
leading  edge  of  solute  and  increase  simulated  concentrations 
at  the  observation  point. 

All  models  for  which  DDA  output  was  available  con¬ 
sistently  underpredicted  the  observed  concentrations  (fig.  8). 
CALF  output  for  DDA  was  unstable  and  included  negative 
values;  however,  prior  simulations  with  CALF  accurately 
simulated  atrazine,  DEA,  and  DDA  at  the  White  River  Basin 
site  (Bayless,  2001).  The  current  study  used  values  of  moisture 
content  at  field  capacity  and  wilting  point  from  an  empirically 
derived  water-retention  curve,  to  provide  consistent  estimates 
for  models  that  use  different  water-retention  functions.  These 
values  differed  from  values  used  previously  in  the  White  River 
Basin  and  might  have  adversely  affected  CALF  predictions. 
Testing  of  CALF  and  LEACHP  included  tracking  a  second 
daughter  product  of  atrazine  (HYA),  but  the  predicted  concen- 
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Figure  6.  Atrazine  concentration  with  time  at  White  River  Basin, 
Indiana.  Atrazine  applied  on  April  28, 1994  (Julian  day  118).  Units 
are  micrograms  per  liter  of  water. 


Figure  7.  Desethylatrazine  concentration  with  time  at  White 
River  Basin,  Indiana.  Atrazine  applied  on  April  28, 1994  (Julian  day 
1 1 8).  Units  are  micrograms  per  liter  of  water. 
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Figure  8.  Didealkylatrazine  concentration  with  time  at  White 
River  Basin,  Indiana.  Atrazine  applied  on  April  28, 1994  (Julian  day 
1 1 8).  Units  are  micrograms  per  liter  of  water. 


trations  are  consistently  less  than  0.005  i.ig/L;  therefore,  HYA 
results  are  not  plotted. 


Model  Summaries 

The  five  modelers  conducting  these  simulations  have 
provided  a  brief  summary  of  each  model,  based  on  their  expe¬ 
riences  with  the  two  data  sets.  The  summaries  emphasize  gen¬ 
eral  model  capabilities  and  key  strengths  and  weaknesses  that 
emerged  during  these  evaluations.  Advantages  and  disadvan¬ 
tages  of  each  model  are  discussed  first  for  the  complex  models 
(LEACHP,  HYDRUS2D,  RZWQM,  and  VS2DT),  and  then  for 
the  simple  models  (CALF,  GLEAMS,  and  PRZM). 


Complex  models 

LEACHP 

Advantages 

The  Leaching  Estimation  and  Chemistry  Model 
(LEACHM)  is  a  DOS-based,  modular  model  for  independently 
simulating  transport  of  heat,  inorganic  chemicals,  pesticides, 
and  nutrients  in  the  unsaturated  zone.  The  LEACHP  (Leach¬ 
ing  Estimation  and  Chemistry  Model  for  Pesticides)  module 
of  LEACHM  simulates  pesticide  transport  using  Richards’ 
equation  for  water  flow  in  the  unsaturated  zone.  Input  data 
are  entered  into  a  single  file  that  can  be  edited  with  any  line 
editor  (vi.  Notepad,  etc.).  The  input  file  is  logically  arranged 
and  includes  extensive  comments;  users  can  see  all  data 
entered  for  a  single  simulation  in  one  long  page,  as  opposed  to 
the  multiple  files,  windows,  and  menus  that  were  used  in  the 
other  codes  tested.  LEACHP  has  no  databases  to  supply  input 


information;  however,  demonstration  input  files  are  supplied 
for  each  module.  LEACHP  can  estimate  characteristic  curve 
and  unsaturated  conductivity  functions  from  grain-size  data  or 
can  utilize  user-supplied  parameters  for  the  Campbell  equa¬ 
tion.  LEACHP  can  simulate  up  to  eight  pesticides  separately 
or  in  degrading  succession.  The  LEACHM/LEACHP  code 
and  documentation  are  available  for  free  at  http.V/www.scieng. 
flinders.edu.au/cpes/people/hutson _j/leachweb.html  (accessed 
on  December  16,  2003). 

Disadvantages 

The  modularity  (and  simplicity)  of  the  LEACHM  code 
makes  it  user  friendly,  but  also  results  in  exclusion  of  some 
processes  and  environmental  data  that  have  bearing  on  pesti¬ 
cide  fate  and  transport  simulations.  For  example,  the  capacity 
version  of  LEACHM  divides  water  into  mobile  and  immobile 
zones,  but  it  is  not  a  true  dual-porosity  model  because  the  two 
zones  cannot  be  independently  controlled;  and,  the  model 
excludes  pesticides  and  soils  databases.  LEACHM  pays  lim¬ 
ited  attention  to  farm-management  practices  and  crop  develop¬ 
ment.  Simulations  of  nutrient  and  organic-residue  distribution 
and  inorganic  soil  and  water  chemistries  can  be  examined  in 
other  LEACHM  modules;  however,  because  those  modules  are 
uncoupled  from  the  LEACHP  module,  their  influence  on  pesti¬ 
cide  fate  cannot  be  directly  simulated.  Similarly,  water  flow 
in  LEACHP  is  simulated  entirely  by  solution  of  the  Richards’ 
equation;  a  macropore  simulating  module,  CALF,  is  part  of 
the  LEACHM  code  but  executes  as  a  stand-alone  module  (see 
below).  There  is  no  graphical  output  of  LEACHP  results  and 
this  slows  down  manual  calibration.  LEACHP  does  not  have 
built-in  capability  to  perform  parameter  estimation  or  uncer¬ 
tainty  analysis. 


HYDRUS2D 


Advantages 

HYDRUS2D  has  extensive  capabilities  for  process-ori¬ 
ented  modeling.  Because  the  model  is  fairly  extensive,  it  may 
require  a  bit  more  time  to  learn  than  other  models.  Available 
documentation  consists  of  help  within  the  graphical  user 
interface  (GUI),  and  a  manual  (available  online).  The  help 
files  for  the  program  are  indexed  and  searchable,  but  only 
cover  the  basics  of  the  model  GUI  and  parameter  set.  More 
extensive  information  is  available  within  the  manual,  which 
is  not  searchable.  User  forums  and  a  list  of  frequently  asked 
questions  are  available  online  at  ]/ittp://www.pc-progress.\ 


cz/Fr_Hydrus. htrrj  (accessed  on  January  22,  2004). 

There  are  two  versions  of  HYDRUS:  HYDRUS  ID  and 
HYDRUS2D.  The  one-dimensional  (1-D)  version  is  available 
online  for  free.  The  model  simulates  nearly  identical  processes 
as  the  2-D  version,  with  some  added  features.  The  solution 
time  is  faster  for  1  -D  problems  and  the  developers  suggest 
that  the  1  -D  version  is  more  appropriate  for  certain  types  of 
problems  than  the  2-D  version.  The  2-D  version  is  available 
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for  nominal  cost,  but  the  optional  grid  generator,  which  is 
important  for  complex  2-D  problems,  does  cost  extra.  The 
site  licensing  policy,  however,  is  that  three  individual  licenses 
at  one  site  constitute  a  site  license,  and  additional  copies  of 
the  model  are  free.  The  GUI’s  of  the  two  versions  are  nearly 
identical,  and  if  a  user  knows  one  version,  the  other  will  be 
easy  to  use. 

HYDRUS2D  has  many  features  that  may  not  be  available 
in  other  Richards’  equation-based  models.  For  example,  there 
are  options  to  automatically  calculate  a  tortuosity  factor  for 
the  diffusion  coefficient,  an  option  to  simulate  hysteresis  in  the 
retention  curves,  and  options  to  include  physical,  non-equi¬ 
librium  in  solute  transport  (mobile-immobile  zones).  Addi¬ 
tionally,  HYDRUS2D  has  coupled  heat  and  water  transport 
capability.  These  features  may  be  helpful  for  testing  concep¬ 
tual  models  of  solute  transport. 

HYDRUS2D  has  built-in  soil  catalogs  with  default 
parameters  for  different  soil  types.  There  is  also  an  option 
of  using  neural  network  predictions  to  get  Rosetta  estimates 
of  water  retention  parameters  and  hydraulic  conductivity. 
HYDRUS2D  has  a  database  of  suggested  values  for  root  water 
uptake  for  different  plants,  and  most  of  the  common  crops  are 
included  in  this  database. 

The  output  graphics  are  clear  and  include  2-D,  X-Y 
graphs,  contour  maps,  and  spectral  maps.  These  can  appear 
on-screen  or  be  exported  as  images.  It  is  also  possible  to  get 
graphs  of  any  variable  along  an  arbitrary  line  across  the  model 
domain. 

Additional  strengths  of  the  model  include:  convenient 
input  for  precipitation  and  ET  on  a  daily  basis;  a  user-friendly 
GUI;  tabular  outputs  of  retention  and  conductivity  curves; 
calculations  of  mass  budgets  in  arbitrary  regions  at  any  depth; 
flexible  and  robust  characteristic  curve  models;  ability  to 
simulate  five  solutes  either  independently  or  in  a  reaction 
chain;  and  built-in  inverse  modeling  capability. 

Disadvantages 

HYDRUS2D  has  several  limitations.  The  assignment  of 
boundary  conditions  is  a  little  complicated  and  confusing.  For 
a  given  node,  the  assigned  boundary  type  cannot  change  dur¬ 
ing  a  single  simulation.  For  example,  an  atmospheric  boundary 
node  (which  is  able  to  switch  between  ET  and  precipitation) 
cannot  change  to  a  seepage  face. 

When  HYDRUS2D  does  not  converge,  it  often  returns  a 
floating  point  error  with  no  explanation.  Resolution  of  prob¬ 
lems  with  convergence  was  sometimes  time  consuming  and 
awkward. 

Time-step  control  parameters  in  the  model  are  set  for  the 
entire  simulation  and  not  for  each  specified  variable  boundary 
condition  period.  Sudden  changes  in  conditions  (e.g.,  heavy 
rain  following  a  prolonged  dry  period)  may  require  manipula¬ 
tion  of  the  boundary  conditions  (e.g.,  spreading  out  the  rain 
over  a  longer  period)  to  get  the  model  to  converge. 

HYDRUS2D  can  track  five  solutes,  but  the  degradation 
pathways  must  be  in  linear  degradation  sequences.  Thus,  the 


model  cannot  simulate  more  than  one  daughter  of  the  parent 
compound  in  the  same  simulation. 

Finally,  the  root  water  uptake  information  is  based  on  a 
single  crop-growing  season.  If  the  simulation  period  includes 
more  than  one  crop  cycle,  the  model  needs  to  be  run  as  succes¬ 
sive  simulations.  The  model  does  not  allow  time- variable  root 
spatial  distribution. 

RZWQM 

Advantages 

Root  Zone  Water  Quality  Model  (RZWQM)  is  a  multi¬ 
process,  Windows-compatible  program  for  simulating  heat, 
nutrient  and  pesticide  transport  in  the  unsaturated  zone  under 
agricultural  influence.  Soil  heat  flux  during  water  redistribu¬ 
tion  is  calculated  with  an  equation  that  includes  terms  for 
diffusive  transport,  convective  heat  transport  by  liquid,  and 
latent  heat  transfer  by  vapor.  Heat  advection  during  infiltra¬ 
tion  events  is  simulated  with  a  separate  heat  balance  equa¬ 
tion.  Three  input  files  for  climate  data  must  be  created  prior 
to  model  execution.  The  climate  files  can  be  auto-generated, 
based  on  probabilistic  analysis  of  long-term  historical  records 
from  meteorological  stations,  or  the  user  can  create  them  from 
on-site  measurements.  All  other  information  relevant  to  the 
simulations  is  entered  through  drop-down  menus.  Pesticide- 
property,  farm-management,  and  crop-development  databases 
can  provide  default  information  where  measured  data  are 
unavailable. 

The  model  is  very  flexible  with  the  ability  to  simulate 
more  agricultural  processes  than  any  other  code  tested.  A  fully 
integrated,  physically  based  macropore-simulation  capability 
is  a  fairly  unique  property  of  RZWQM.  The  model  has  the 
ability  to  simultaneously  simulate  transport  of  three  pesti¬ 
cides,  which  may  include  a  parent-daughter-granddaughter 
succession.  RZWQM  contains  modules  to  simulate  overland 
flow,  sediment  routing,  tile  drainage,  inorganic  water  and  soil 
chemistry  (which  can  effect  pesticide  properties),  and  micro¬ 
bial  activity  stimulated  by  organic  residue. 

Model  output  includes  sophisticated  2-D  and  three- 
dimensional  (3-D)  plots  of  time-dependent  input  and  output 
variables  that  can  be  viewed  on  screen  or  saved  in  bitmap  and 
metafile  formats.  Results  are  also  output  in  user-specified 
tabular  files.  The  model  is  well  documented  in  online  indexing 
and  contents  windows,  and  has  an  accompanying  text:  Root 
Zone  Water  Quality  Model  -  Modeling  Management  Effects 
on  Water  Quality  and  Crop  Production  (Ahuja  and  others, 
2000).  The  RZWQM  code  is  available  for  free  from  the  U.S. 
Department  of  Agriculture- Agricultural  Research  Service 
(USDA-ARS)  at  http://gpsr.ars.usda.gov/products/rzwqm.htm 
(accessed  in  January  2004). 

Disadvantages 

The  complexity  of  the  model  results  in  more  parameters 
and,  therefore,  greater  data-input  needs.  A  parameter  estima- 
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tion  and  uncertainty  analysis  package  would  be  a  welcome 
addition.  The  crop-development  database  is  incomplete.  In  the 
tested  version  of  the  code,  the  maximum  simulation  depth  was 
3  m;  but  this  has  since  been  extended  to  30  m  by  the  USDA- 
ARS.  RZWQM  cannot  easily  simulate  macropore  flow  in  dry, 
sandy  soils,  because  the  model  requires  ponding  and  overland 
flow  to  induce  macropore  flow.  Macropore  flow,  however,  can 
be  induced  in  sandy  soils  by  specifying  a  surface  crust  with  a 
low  hydraulic  conductivity.  There  is  no  user-defined  dispersiv- 
ity  in  RZWQM.  Mechanical  dispersion  (function  of  dispersiv- 
ity)  is  simulated  for  mobile  water  by  two-stage,  partial  piston 
displacement,  and  molecular  diffusion  (the  other  component  of 
hydrodynamic  dispersion)  is  explicitly  included  and  parame¬ 
terized  by  the  user  (Ahuja  and  others,  2000).  Alternatively,  the 
scale  relation  with  dispersivity  can  be  simulated  with  numeri¬ 
cal  dispersion.  The  RZWQM  reference  text  (Ahuja  and  others, 
2000)  is  not  free,  but  can  be  purchased  from  Water  Resources 
Publications,  Inc. 

VS2DT 

Advantages 

Variably  Saturated  Two  Dimensional  Transport  (VS2DT) 
is  a  suitable  model  for  simulation  of  unsaturated  flow  and 
transport  under  many  circumstances.  The  GUI  is  relatively 
straightforward  and  easy  to  use.  Basic  help  is  available  from 
within  the  GUI.  Printed  documentation  includes  the  original 
manual  for  VS2D  (flow  only),  the  manual  for  VS2DT  and 
VS2DH  (solute  and  heat  transport,  respectively)  and  a  manual 
for  the  GUI.  All  are  free,  but  only  the  GUI  manual  is  available 
in  electronic  form. 

The  main  advantages  of  VS2DT  are  that  it  is  simple 
to  use  and  easy  to  learn.  The  model  itself  is  free  and  avail¬ 
able  online  from  | http ://water.usgs. gov/softwa re/vs 2 di.li tm i| 
(accessed  on  January  22,  2004).  Several  example  problems 
are  included  with  the  code  to  acquaint  the  user  with  the  GUI. 
When  the  pre-programmed  characteristic  curve  models  are  not 
adequate  for  a  user’s  needs,  any  shape  of  characteristic  curve 
can  be  modeled  by  using  the  tabular  data  input  option.  Another 
feature  of  the  model  is  the  flexibility  of  boundary  conditions 
and  iteration  parameters.  The  assigned  boundary  type  can 
be  changed  for  any  cell  at  any  time.  Likewise,  the  iteration 
parameters  are  assigned  for  each  recharge  period,  so  that  the 
user  can  assign  smaller  time  steps  when  there  is  a  sudden 
change  in  conditions  (e.g.,  rainfall  after  a  prolonged  dry  spell). 
Root  depth  and  activity  can  also  be  varied  with  time.  VS2DT 
has  built-in  databases  with  textural  classes  and  new  classes 
can  be  imported  from  text  files.  The  model  has  22  generic  soil 
types  with  default  parameters  for  these  types. 

The  VS2DT  GUI  has  runtime  graphics  to  show  various 
output  variables,  and  the  graphs  at  any  instant  can  be  exported 
and  saved  as  bitmaps.  Also,  text  files  with  the  output  data  can 
be  imported  into  graphics  packages  such  as  Excel. 


Disadvantages 

The  model  tracks  a  single  compound  and  cannot  simu¬ 
late  the  formation,  transport,  and  fate  of  degradates.  Prefer¬ 
ential  flow  can  be  simulated  by  adjusting  Ks  t  or  by  using  a 
bimodal  retention  function  entered  via  tabular  data,  but  there 
is  no  option  to  simulate  storage  in  immobile  zones.  There 
is  no  built-in  option  to  simulate  hysteresis.  As  compared  to 
HYDRUS,  VS2DT  has  fewer  features  such  as  built-in  data¬ 
bases,  parameter  estimation  routines,  and  options  for  graphical 
display  of  results. 

The  current  version  of  VS2DT  has  no  cut  and  paste 
option  for  the  GUI,  and  all  entries  must  be  entered  manu¬ 
ally.  For  simple  scenarios  this  may  not  be  an  issue,  but  for  a 
scenario  with  300  recharge  periods,  the  GUI  is  impractical. 
Typically,  a  user  can  use  the  GUI  for  much  of  the  problem, 
but  then  must  edit  the  text  input  fde  to  enter  a  lengthy  list  of 
recharge  period  information. 

VS2DT  cannot  have  precipitation  and  evaporation  occur¬ 
ring  simultaneously  in  the  same  recharge  period.  Average  rain¬ 
fall  intensity  is  needed  to  predict  the  duration  of  each  rainfall 
event  and  following  evapotranspiration  event.  Although  the 
quality  of  output  graphics  is  good,  the  scale  does  not  automati¬ 
cally  adjust  to  accommodate  changing  values  such  as  pressure 
and  concentration  during  a  simulation. 

Simple  models 

CALF 

Advantages 

The  CALculates  Flow,  leaching  of  pesticides  in  field 
soils  (CALF)  model  is  an  option  within  LEACHM  (the  family 
of  models  that  includes  LEACHP)  that  simulates  water  flow 
using  a  tipping-bucket  approach,  and  is  intended  to  represent 
macropore  transport  of  pesticides.  Input  data  are  entered  into  a 
single  file  that  can  be  edited  with  any  line  editor  (vi,  Notepad, 
etc.).  Similar  to  LEACHR  the  strength  of  CALF  is  the  simplic¬ 
ity  of  the  data  input  file,  which  is  logically  arranged  and  well- 
commented;  users  can  see  all  data  entered  for  a  simulation 
in  a  single  file.  CALF  contains  no  databases  to  supply  input 
information,  however,  a  sample  input  file  is  provided  with  the 
code  and  documentation.  CALF  can  examine  the  transport  of 
up  to  eight  pesticides  or  pesticide  transformation  compounds. 
A  single  toggle  switches  from  CALF  to  LEACHP,  making  it 
easy  to  examine  transport  by  the  two  water- flux  algorithms. 

Disadvantages 

CALF  is  a  simple  model  in  concept  and  was  able  to 
simulate  atrazine,  DEA,  and  HYA  in  the  current  study.  CALF 
output  for  DDA  was  unstable  and  included  negative  values. 
This  was  not  considered  overly  limiting,  however,  because 
CALF  accurately  simulated  atrazine,  DEA,  and  DDA  at  the 
same  site  in  an  earlier  study  (Bayless,  2001). 
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The  model  is  limited  to  linear  sorption  and  is  strongly 
affected  by  the  user-specified  field  capacity.  Like  LEACHP, 
there  are  no  graphic  outputs  and  plotting  data  must  be 
extracted  from  tabular  files.  The  uncoupled  solution  of  Rich¬ 
ards’  equation  (LEACHP)  for  flow  in  the  soil  matrix  and  the 
Addiscott  option  for  preferential  flow  makes  interpretation  of 
results  difficult  because,  under  field  conditions,  these  pro¬ 
cesses  can  occur  simultaneously.  CALF  does  not  have  built-in 
code  to  perform  parameter  estimation  or  uncertainty  analysis. 

GLEAMS 

Advantages 

Groundwater  Loading  Effects  of  Agricultural  Manage¬ 
ment  Systems  (GLEAMS)  is  a  DOS-based  continuous  simula¬ 
tion,  field-scale,  modular  model  based  on  the  Chemicals, 
Runoff  and  Erosion  from  Agricultural  Management  Systems 
(CREAMS)  model.  GLEAMS  is  incorporated  into  USDA’s 
Soil  and  Water  Assessment  Tool  (SWAT)  watershed  model, 
which  facilitates  simulation  of  UZ  processes  to  watershed 
scales.  GLEAMS  is  relatively  simple  to  learn,  and  simplified 
DOS-based  GUI  files  ease  the  building  of  text  input  param¬ 
eter  files  for  the  hydrology,  erosion,  nutrient,  and  pesticide 
modules.  The  erosion  module  is  comprehensive,  allowing  the 
detailed  input  of  the  topographic  features  of  an  area  for  use 
in  predicting  sediment  yield.  A  total  of  366  pesticides  can  be 
simulated  simultaneously  and  a  pesticide  parameter  database 
provides  information  for  336  possible  chemicals.  Application 
running  time  is  not  a  concern,  averaging  a  few  seconds  per 
simulation  and  running  efficiently  with  few  unnatural  termina¬ 
tions  during  simulations.  Output  is  available  for  a  number  of 
different  daily,  monthly,  and  annual  variables  in  both  standard 
and  metric  units,  for  both  calibration  and  display  purposes. 

Disadvantages 

Climate  data  files  for  daily  temperature  and  precipitation 
require  a  specific  format  that  is  tedious  to  produce  for  model 
simulations.  All  other  climate  data  are  entered  as  monthly 
averages.  Some  parameters,  including  saturated  hydraulic  con¬ 
ductivity  and  porosity,  are  limited  to  intervals  specified  in  the 
input  module  GUI.  Simulation  conditions  are  limited  to  five 
soil  horizons  and  confined  to  the  root  zone,  which  is  limited 
to  a  maximum  depth  of  1.52  m.  Irrigation  estimates  can  be 
calculated  using  an  internal  process  that  balances  soil-moisture 
values,  but  will  not  allow  for  the  input  of  specific  irrigation 
values.  The  inclusion  of  these  irrigation  totals  within  precipita¬ 
tion  values  allows  for  a  possible  spike  in  runoff  that  would  not 
occur  during  normal  irrigation  practice.  While  the  pesticide 
and  nutrient  modules  of  the  application  can  be  excluded  from 
the  simulation,  the  erosion  module  is  mandatory  for  each  sim¬ 
ulation,  which  may  be  beyond  the  desired  scope  of  the  project. 
Simulation  output  is  available  for  many  different  variables, 
but  some  common  variables  such  as  soil  moisture  and  matric 
pressure  are  absent. 


PRZM 

Advantages 

The  Pesticide  Root  Zone  Model  (PRZM)  is  capable  of 
simulating  complex  scenarios  such  as  the  metabolic  activity  of 
microbial  populations,  pesticide  application  techniques,  ero¬ 
sion,  runoff,  plant  uptake,  and  agricultural-management  prac¬ 
tices.  PRZM  also  simulates  pesticide  decay  in  the  dissolved, 
adsorbed,  and  vapor  phases.  There  is  no  maximum  limit  for 
depth  of  simulation,  and  the  ability  to  vary  horizon  thickness 
is  useful  for  representing  detailed  soils  data.  PRZM  produces 
time-series  data  of  fluxes  at  a  given  point  in  the  soil  profile  in 
columnar,  ASCII  format  that  can  be  processed  fairly  quickly. 

Disadvantages 

Model  output  is  extremely  sensitive  to  horizon  (grid  cell) 
thickness,  which  is  arbitrarily  set  by  the  model  user.  Model 
documentation  provides  minimal  guidance  on  selection  of  the 
number  and  thicknesses  of  horizons. 

The  user’s  manual  is  sometimes  difficult  to  follow  and 
focuses  more  on  theory  than  on  practical  guidance  in  using 
the  model.  Although  there  is  extensive  discussion  of  theory,  it 
is  sometimes  unclear  how  or  whether  the  model  simulates  the 
processes  being  described.  PRZM  was  written  in  FORTRAN 
in  the  early  1980s  and  uses  formatted  ASCII  files  for  input. 

As  a  consequence,  learning  how  to  set  up  and  operate  PRZM 
can  be  difficult.  For  example,  learning  how  to  format  the  main 
input  file  can  be  time  consuming. 


Summary 

Seven  unsaturated-zone  contaminant-transport  models 
were  compared  for  ease  of  use  and  model  capabilities.  The 
models  were  tested  with  two  data  sets  in  cold  simulation 
mode.  Proper  calibration  was  beyond  the  scope  of  the  study, 
but  differences  in  model  output  provided  insights  into  impor¬ 
tant  differences  in  model  capabilities. 

Among  complex  models,  RZWQM,  HYDRUS2D,  and 
VS2DT  were  the  easiest  to  use  because  of  user-friendly 
GUIs  and  extensive  documentation.  Of  these,  RZWQM  and 
HYDRUS2D  can  simulate  several  parent  compounds  and 
degradates.  RZWQM  can  explicitly  simulate  water  and  solute 
flux  in  macropores,  a  potentially  important  feature  for  sites 
with  structured  soils.  Both  HYDRUS2D  and  VS2DT  can 
simulate  water  and  solute  flux  in  two  dimensions,  which  might 
be  needed  at  sites  near  streams  or  drains.  Although  the  tested 
version  of  RZWQM  has  a  maximum  simulation  depth  of  3 
m,  the  USDA-ARS  has  since  extended  this  depth  to  30  m. 
VS2DT  simulates  the  transport  and  fate  of  a  single  compound. 
Although  LEACHP  has  no  GUI,  it  is  logically  arranged  and 
relatively  easy  to  use.  The  Richards’  equation  version,  how¬ 
ever,  lacks  preferential  flow  capability. 
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Considering  all  models,  only  HYDRUS2D  predicted 
both  surface-retention  of  bromide  and  the  deeper  peak  at  the 
Merced,  CA  site,  based  on  a  mobile-immobile  water-flow 
simulation.  RZWQM  cannot  easily  simulate  macropore  flow 
in  dry,  sandy  soils,  because  the  model  requires  ponding  and 
overland  flow  to  induce  macropore  flow.  RZWQM  predictions 
of  atrazine  (White  River  Basin)  were  closest  to  the  measured 
values.  RZWQM  simulated  macropore  flow  at  this  site,  which 
is  wetter  than  the  Merced  site  and  has  heavily  structured  soils 
with  visible  surface  cracks. 

Simple  models  can  be  used  for  scale-up  and  mass 
balances  at  the  stream-basin  scale.  Among  simple  models, 
CALF  (the  tipping-bucket  version  of  LEACHP)  is  simple  in 
concept  and  easy  to  use,  and  has  preferential  flow  capability. 
GLEAMS  is  easy  to  use,  simulates  effects  of  agricultural  prac¬ 
tices,  and  is  incorporated  into  USDA’s  Soil  and  Water  Assess¬ 
ment  Tool  (SWAT)  watershed  model.  GLEAMS,  however,  has 
a  maximum  simulation  depth  of  1.5  m,  and  the  authors  are 
unaware  of  any  plans  to  extend  it.  The  depth  limitation  creates 
a  potential  disconnect,  for  mass  balance  purposes,  between 
unsaturated-zone  transport  and  the  water  table.  PRZM  has 
no  depth  limitation  and  can  simulate  a  variety  of  processes, 
including  microbial  population  growth,  plant  uptake,  and 
agricultural-management  practices.  PRZM,  however,  can  be 
difficult  to  use  and  the  documentation  is  hard  to  follow.  Addi¬ 
tionally,  PRZM  is  sensitive  to  numerical  grid-cell  thickness, 
which  determines  effective  dispersivity. 

Among  tipping-bucket  models,  CALF  predicted  the 
highest  values  of  atrazine  and  DEA  at  the  measurement  depth 
of  1.5  m  at  the  White  River  Basin  site.  CALF  includes  the 
Addiscott  flow  option  for  preferential  flow,  and  also  accepts 
user-specified  dispersivity.  GLEAMS  and  PRZM  predicted 
negligible  concentrations  of  atrazine  and  degradates  at  a  depth 
of  1.5  m.  Low  pesticide  concentrations  predicted  by  PRZM 
apparently  resulted  from  the  numerical  layering  scheme  (1 
mm  in  the  upper  10  cm  of  soil  and  1  cm  at  greater  depths). 
Increasing  dispersion  (via  grid  size  in  this  case)  would  acceler¬ 
ate  movement  of  the  leading  edge  of  solute  and  would  increase 
simulated  concentrations  at  the  observation  point. 

The  unsaturated-zone  model  evaluation  provided  insights 
into  model  capability  and  potential  shortcomings.  Capabili¬ 
ties  assessed  include  the  number  of  dimensions,  macropore 
flow,  multiple  compounds,  and  GUIs.  For  some  models,  clear 
limitations  in  simulation  depth,  effects  of  grid  cell  size,  and 
GUIs  for  input  and  output  were  found.  In  general,  the  study 
permitted  evaluation  of  these  factors  in  the  context  of  antici¬ 
pated  modeling  needs  for  the  unsaturated-zone  component  of 
NAWQA  ACT  studies. 
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